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I. OVERVIEW 

An exact simulation of quantum systems on a classical computer is computationally hard. 
The problem lies in the dimensionality of the Hilbert space needed for the description of 
our system that in fact grows exponentially with the size of this system. No matter if we 
simulate the dynamics or calculate some static property e.g. the energy, this limitation is 
always present. Richard Feynman came up with an alternative to the classical simulation 
[1]. His idea was to convert the aforementioned drawback of quantum systems into their 
benefit. He suggested to map the Hilbert space of a studied quantum system on another one 
(both of them being exponentially large) and thus to efficiently simulate a quantum system 
on another one (i.e. on a quantum computer). 

This was the original idea of quantum computers. There is no doubt that quantum com- 
puting is nowadays a well-established discipline of computer science. Apart from the efficient 
simulation of quantum systems [2HS1 , other interesting problems where quantum computers 
could beat their classical counterparts have been discovered. The most famous examples are 
integer factorization for which quantum computers supply an exponential speedup [71 18] or 
database search [9]. However, for our purposes the efficient (polynomially scaling) quan- 
tum algorithm of Abrams and Lloyd for obtaining eigenvalues of local Hamiltonians [10] is 
particularly important. 

The first paper connecting quantum computation and chemistry was published by Lidar 
and Wang [TT] and concerned the efficient calculations of thermal rate constants of chem- 
ical reactions. This work in fact founded the new field of computational chemistry that 
is the subject of this book, namely the "computational chemistry on quantum computers" . 
Aspuru-Guzik et al. in their seminal article [12] reduced the number of quantum bits (qubits) 
needed by the Abrams and Lloyd's algorithm [Tfl] and applied it to molecular ground state 
energy calculations. Since these two pioneering works, other papers involving energy cal- 
culations of excited states [13J, quantum chemical dynamics [H], calculations of molecular 
properties and geometry optimizations [15], state preparations [161 IH] or global minima 
search [T8| were published. The list of all chemical applications for quantum computers is 
quite rich and an interested reader can find a comprehensive review in |19j . 

Aspuru-Guzik et al. [12] also proposed that quantum computers with tens of (noise 
free) qubits would already exceed the limits of classical full configuration interaction (FCI) 
calculations. This is in contrast to other quantum algorithms, e.g. the Shor's algorithm 
[3 IE] for integer factorization would for practical tasks in cryptography require thousands 
of qubits. For this reason, calculations and simulations of quantum systems will belong to 
the first practical applications of quantum computers. Recent proof-of-principle few-qubit 
experiments covering energy calculations of the hydrogen molecule [201 [21] or Heisenberg spin 
model [22] and the simulation of a chemical reaction dynamics [23] confirm that interesting 
applications might be just behind the door. 

The aim of this article is to review quantum algorithms for exact molecular energy calcula- 
tions (within a finite one-particle basis set), i.e. quantum analogues of the FCI calculations. 
On a classical computer, a computational cost of the FCI method scales exponentially with 
the size of the system. This fact stems from the dimension of the Hilbert space in which we 
diagonalize the Hamiltonian matrix and it is the reason why this method is limited only to 
the smallest systems (diatomics, triatomics). For example, in the non-relativistic case, the 
number of Slater determinants that build up the FCI wave function for a closed-shell system 



with n electrons in m orbitals is equal to 

A^non-rel. = (^ J/ 2 J " (l) 

It is evident that this number grows into huge values with increasing m very quickly. On a 
quantum computer on the other hand, it has been shown |20l El] and will be discussed later 
in the article that the FCl cost has a polynomial scaling [(9(m^)], therefore it is exponentially 
faster. 

The structure of this article is as follows: First we very briefly give the necessary basics 
of quantum computing including the quantum Fourier transform (QFT), the phase estima- 
tion algorithm (PEA) and their semiclassical variants (Section O. Despite being textbook 
topics [25j, we believe that it may be convenient for readers from the quantum chemistry 



community to make this article more self-contained. Section III presents details of the 



quantum algorithm for the FCI method (qFCI). In Section IV, we show its applications to 
non-relativistic computations, namely the computations of the ground and excited state en- 
ergies of the methylene molecule (CH2). We then generalize this approach to the relativistic 
4-component (4c) calculations and apply it to the problem of spin-orbit splitting in the SbH 
molecule in Section |Vl 



II. QUANTUM COMPUTING BACKGROUND 

In this section we briefly mention the necessary quantum computing background. For 
a more detailed description, we refer the reader e.g. to the excellent book by Nielsen and 
Chuang [2^. 

A. Quantum Fourier transform 

The classical discrete Fourier transform takes as an input a vector of complex numbers 
(xo, . . . , xn-i) and outputs the elements of another vector {yo, . . . , yN-i) according to the 
equation 
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Similarly, the quantum Fourier transform (QFT) operates on an orthonormal basis of n 
qubits: |0) . . . 12"" — 1) and is defined as an operator J/qft 

N-l 

UQFTlk) = -7^J2 e'^'^'^'^lj)^ N = 2\ (3) 

i=o 

where the kets are numbered by a binary representation of integers 

|j) = bn---Jl), J»G{0,1}, 

i=l 
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Figure 1: The quantum Fourier transform circuit. Note that qubits of the result are in a 
reversed order after the apphcation of this circuit. 



It can be shown that the QFT is an unitary operator [25j. 

The advantage of the QFT is that it can be performed with just C(n^) operations (quan- 
tum gates). This is in sharp contrast to the classical fast Fourier transform (FFT) with 
the scaling 0{N\ogN = n2"'). Without the derivation, which can be found e.g. in ^25] . 
the quantum circuit for the QFT is shown in Figure [I] H (Hadamard) and Rj gates are 
represented by the following matrices 



H 



1 

V2 



Rj 



1 

e2-/2^ 



(5) 



It should be noted that even though the QFT can be done exponentially faster than the 
FFT, it cannot be used as an efficient straightforward replacement of the Fourier transform 
itself. It would indeed require to prepare an arbitrary state of n qubits and also measure 
all of the complex amplitudes at the end, which cannot be done efficiently. Nevertheless 
the QFT is a key part of the phase estimation algorithm [25] (contained also in the Shor's 
factoring algorithm [TJ E]) as will be shown in Section II C 



B. Semiclassical approach to quantum Fourier transform 

The QFT circuit from Figure [l] can in fact be greatly simplified using the semiclassical 
(measurement based) approach [26]. Notice that gates acting on each qubit [except the first 
(top most) and the last ones where the corresponding parts are missing] obey the general 
structure: first, Rj gates controlled by previous qubits are applied, then the Hadamard gate 
is applied, and finally they serve as control qubits for the subsequent ones. Because the 
state of each qubit does not change after the application of the Hadamard gate, we can in 
fact do the measurement just after this gate. Instead of employing controlled Rj gates, we 
then apply only the corresponding one qubit gates depending on the results of individual 
measurements. Moreover, all Rj gates acting on a kth qubit can be merged into a single 
rotation gate 
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whose angle ujk depends on the results of previously measured qubits (gj) according to the 



formula 
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Figure 2: Simplified, measurement based circuit for the A;th qubit of the QFT. 



Figure |2] shows the semiclassical QFT circuit pattern which is the same for all qubits. 
The big advantage of the aforementioned approach is that we have actually replaced two 
qubit gates by one qubit ones (controlled by a classical signal). This technique is especially 
useful in connection with the phase estimation algorithm where it leads to the formulation 



of its iterative version (IPEA, Section II D) 



C. Phase estimation algorithm 

The phase estimation algorithm (PEA) [25] is a quantum algorithm for obtaining the 
eigenvalue of an unitary operator f/, based on a given initial guess of the corresponding 
eigenvector. Since an unitary f/ can be written as f/ = e*"^, with H Hermitian, the PEA 
can be viewed as a quantum substitute of the classical diagonalization. 

Suppose that \u) is an eigenvector of f/ and that it holds 



lJ\u) 



^2-Ki4>\ 



m), 



e(0,i), 



(8) 



where is the phase which is estimated by the algorithm. Quantum register is divided into 
two parts. The first one, called the read-out part, is composed of m qubits on which the 
binary representation of is measured at the end and which is initialized to the state lO)®™. 
The second part contains the corresponding eigenvector \u) . 
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Figure 3: The PEA circuit with the highlighted part corresponding to the inverse QFT. 
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The PEA quantum circuit is shown in Figure [3j The apphcation of Hadamard gates on 
all qubits in the first part of the register gives 



2'"-! 

Next, after the application of a sequence of controlled powers of U, the register is transformed 
into 



g) = i E u^\^)\^) = 4^ E ^''''''\j)\^)- (10) 



o-m f J /Or 

The heart of the PEA is the inverse quantum Fourier transform (QFT''', highlighted in Figure 
[3]) performed on the read-out part of the register which is transformed to |2™0)|'u). 
If the phase can be expressed exactly in m bits 

= 0.0102 . . . 0^ = ^ + ^ + . . . + ^, 0. G {0, 1}, (11) 

it (and consequently the eigenvalue) is recovered with unity probability by a measurement 
on the first part of the quantum register. 

The situation is more complicated when cannot be expressed exactly in m bits. Then 
we can write 

= + (52-'", (12) 

where = 0i02 . . . 0m denotes the first m bits of the binary expansion and b : < (5 < 1 is 
a remainder. The closest m-bit estimators of correspond to either (rounding down) or 
0+2-™ (rounding up). When we label the probabilities of measuring these two estimators by 
-Pdown and Pup, it can be shown (e.g. [2Z1) that the sum Pdown + -Pup decreases monotonically 
with increasing m. The dependence of Pdown and Pup on (5 for m = 20 is presented in Figure 
El In the limit m — ?■ 00, the lower bound reads [27] 

Pdown(<5 = 1/2) + Pup(5 = 1/2) = 1 + 1 > 0.81. (13) 

If the desired eigenvector is not known explicitly (as is typically the case in quantum 
chemistry), we can start the algorithm with an arbitrary initial guess vector I?/)), which can 
be expanded in terms of eigenvectors of f/ 



J^Cilwi). (14) 



The probability of obtaining the exact m-bit 0j is due to linearity of the algorithm |cjp. 
It is important to note that the initial guess does not influence the accuracy of the phase, 
only the probability with which the phase of a particular eigenstate is measured. When 0j 
cannot be expressed in m bits, the lower bound for Pdown + Pup corresponding to 0j is equal 
to 0.81- Icf. 
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Figure 4: The dependence of success probabilities of the PEA on (5 for m = 20. Pdown and 
Pup denote the success probabihties corresponding to rounding the exact phase up/down to 

m binary digits. 



D. Iterative phase estimation algorithm 



Using the semi classical QFT [2S] (Section |IIB ), the PEA circuit can be simplified, having 
only one ancillary qubit in the read-out part of the quantum register. The algorithm then 
proceeds in an iterative manner [iterative phase estimation algorithm (IPEA)]. The /c-th 
iteration of this scheme is presented in Figure |5| Note that as the PEA uses the inverse 
QFT, the angle oj^ (JTl) must be negative now. The algorithm is iterated backward from 
the least significant bits of to the most significant ones, for k going from m to 1. For 
our purposes, the presented IPEA, which is the unitary matrix eigenvalue algorithm, is 
completely adequate, but we would like to note that Wang et al. recently presented a 
modified version of the IPEA capable of finding eigenvalues of non-unitary matrices [28] . 
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Figure 5: The fc-th iteration of the IPEA. The feedback angle depends on the previously 
measured bits, k is iterated backwards from m to 1. 



The IPEA is in fact completely equivalent to the original (multiqubit) PEA [27j. It 
exhibits the same behaviour - decreasing of the success probability when the phase cannot 
be expressed exactly in a particular number of bits. One possibility of a success probability 
amplification is performing more iteration steps (more than is the desired accuracy of 0): 
when extracting m' = m + log(2 + l/2e) bits, the phase is accurate to m binary digits with 
probability at least 1 — e [23] . This method is however not very useful, since implementing 
for large k is the algorithm's bottleneck in a realistic noisy environment j29j . 



U^' 



Another alternative [29] is to repeat the measurement for the least important bits of 
the phase binary expansion. Using the majority voting (for bit value or 1), the effective 
error probability decreases exponentially with the number of repetitions according to the 
binomial distribution. This measurement repetition only for the few least important bits of 
is unfortunately possible only if the exact eigenstates of U are available. 
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a) Maintaining the second part of the quantum register during all iterations (version A) . 
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b) Repeated initial state preparation in each iteration (version B). 
Figure 6: Comparison of the two versions of IPEA, ISP denotes the initial state 



preparation. 

When working with general initial states (as in a practical application to quantum chem- 
istry), two scenarios are possible [30], as shown in Figure |6} Maintaining the second part 
of the quantum register during all iterations and amplification of the success probabilities 
by repeating the whole process number of times is the first possibility. This version was 
denoted as A version of IPEA. The biggest advantage of this approach is that one always 
ends up with one of the eigenstates of f/ in the second part of the quantum register as 
was also the case of the PEA. It happens through successive collapses of the system state 
into the corresponding eigensubspace and is demonstrated on the hydrogen molecule with 
random initial states in Figure [7} The biggest disadvantage that complicates the potential 
physical realization of this scheme is the requirement for a long coherence time of the quan- 
tum register. We would like to note here that when amplifying the success probability by 
repeating the whole process, it must be higher than 0.5 to be sure that we get the energy 
of the right state. This, however, is not necessary for the ground state energy which can 
always be identified by the lowest eigenvalue [221 ED] ■ 

Another possibility is to initialize the second part of the quantum register at every it- 
eration step (B version of IPEA). Every iteration step (not only the least important bits 
of as in [29]) must be repeated and measurement statistics performed. One could other- 
wise possibly mix bits belonging to different eigenvalues in different iterations and obtain 
an unphysical result. The biggest advantage of this approach is avoidance of the long co- 
herence times and therefore potentially easier physical implementation. On the other hand, 
the biggest disadvantage is that no improving of the overlap between the actual state of the 
quantum register and the exact wave function occurs during the iterations and one must 
"fight" the overlap problem at every iteration step. But as will be shown in Section IV on 
the example of the methylene molecule 



[30] , small number of repetitions of each iteration is 
sufficient for amplification of the success probability to unity, when a suitable initial state 
of the quantum register is used. 

At the end of this section, we would like to mention a different way of reducing the 
number of read-out qubits required by the PEA, which was suggested in the seminal work 
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Figure 7: Energies of the four electronic states of H2 in ST0-3G basis which were obtained 
by the qFCI method (IPEA version A) with randomly generated initial guess states. Small 
figure inside presents the increasing overlap between the actual state of the second part of 
the quantum register and the exact wave function for one of random runs of the algorithm 
leading to the ground state. Figure reprinted from [3Dj. 



by Aspuru-Guzik et al. [12]. Their recursive variant of the PEA uses four qubits in the 
read-out part of the quantum register on which the phase (and therefore also the energy) 
is successively improved. It starts with measuring the first four bits of the phase (p. The 
Hamiltonian is then shifted by this reference value and a four-bit estimate of the deviation 
of the phase from the reference one on the half of the interval computed. The procedure is 
iteratively repeated and the overall effect is a gain of one additional bit of at each iteration 
step (thus the same as in the IPEA). In spite of the fact that this method uses four read-out 
qubits instead of one which is used by the IPEA, it is worth mentioning. First of all, it was 
the first iterative version of the PEA applied to the Hamiltonian eigenvalue problem ^12]. 
Secondly, it recovers the energy starting from the most important bits towards the least 
important ones (in contrast to the IPEA), which can be advantageous in certain situations. 



III. QUANTUM FULL CONFIGURATION INTERACTION METHOD 

The PEA/IPEA can in fact be used for efficient computations of eigenvalues of local 
Hamiltonians [TU] . If we take U in the form 



U 



rH 



(15) 



where i^ is a local Hamiltonian and r a suitable parameter assuring being in the interval 
(0, 1), then the algorithm reveals the energy spectrum of H. The whole procedure can be 
simply viewed as a time propagation of a trial wave function followed by the QFT switching 
from the time to energy domain and a measurement projecting out a certain eigenstate. 

In this section, we discuss the application to molecular Born-Oppenheimer electronic 
Hamiltonians. We will start with a mapping of quantum chemical wave functions onto a 
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quantum register (Section [III A). Section IIIB briefly mentions the initial state preparation 



and Section IIIC deals with the crucial part of the algorithm, namely the implementation 
of controlled powers of the exponential of molecular Hamiltonians (controlled "time propa- 
gation" ) . 



A. Mapping of quantum chemical w^ave functions onto quantum register 

Several possible mappings of a quantum chemical wave function onto a register of qubits 
have been proposed. The simplest, but the least economical one in terms of the number 
of employed qubits, is so called direct mapping [12]. In this approach, individual spin 
orbitals are directly assigned to qubits, since each spin orbital can be either occupied or 
unoccupied, corresponding to |1) or |0) states. The inefficiency lies in the fact that it 
actually maps the whole Fock space of the system (states with different number of electrons) 
on the Hilbert space of the quantum register. Relativistic generalization of this approach 
^31j assigns one qubit to one Kramers pair bispinor {A or B, analogy to a and /3 spin in 
non-relativistic treatment). The advantage of the direct mapping stems from the fact, that a 
general factorization scheme, i.e. an algorithm to systematically generate a quantum circuit 



implementing the exponential of a Hamiltonian has been discovered (see Section IIIC 1). 

Compact mappings from a subspace of flxed-electron-number wave functions, spin- 
adapted |12], or even symmetry-adapted wave functions employing the point group [13] 
or double group symmetry [21] to the register of qubits have also been proposed. How- 
ever, general factorization schemes are unfortunately not known for these mappings. The 
factorization to elementary quantum gates can be for small circuits performed either with 
numerical optimization techniques (e.g. with genetic algorithms [32j) or analytically [33] . 
but neither efficiently. Its use is motivated by the need to employ as few qubits as possible 
in today's experimental realizations [201 EI]- 



B. Initial states for the algorithm 

The quantum full conflguration interaction (qFCI) algorithm must be started with some 
initial guess state. Generally, it holds that the closer is the initial guess to the exact wave 
function corresponding to the calculated energy, the higher is the success probability of 
measuring the energy. As was shown in [T71 130], the simplest one-determinantal Hartree- 
Fock guess may not be successful in situations, where correlation (particularly the static 
one) plays an important role. In these situations, initial guess states from more sophisticated 
polynomially scaling methods can be used [e.g. complete active space self consistent fleld 
(CASSCF) method in a limited orbital CAS]. 

Preparing a general initial state (an arbitrarily chosen vector from the Hilbert space 
of n qubits) is a hard task as this vector can contain up to 2" non-zero components in 
general and it cannot be performed efficiently. Fortunately, initial guesses including only 
few determinants in a superposition are sufficient for most purposes of quantum chemistry 
[30] . These states can be prepared e.g. with the procedure described by Ortiz et al. [1] 
which scales as 0{N'^) in the number of determinants A^. Preparation of general molecular- 
like states from the combinatorial space of dimension ( ) corresponding to distributing m 
electrons among n spin orbitals was presented in [T7| . Preparation of many-particle states in 
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a superposition on a lattice which can be then propagated by quantum chemical dynamics 
algorithm [H] was studied in [16]. 



1. Adiabatic state preparation 

A different approach of the initial state preparation is the adiabatic state preparation 
(ASP) method of Aspuru-Guzik et al. [12]. In the ASP method, one slowly varies the 
Hamiltonian of the quantum register, starting with a trivial one and the register in its 
(exactly known) ground state and ending with the final exact one in the following simple 

way 

^ = (1 - S)ifinit + Si^exact S : -^ 1. (16) 

If the change is slow enough (depending on the gap between the ground and the first excited 
state), the register remains in its ground state according to the adiabatic theorem [3^ . 
In the compact mapping, Hinit can be defined to have all matrix elements equal to zero, 
except Hii, which is equal to the (Dirac-)Hartree-Fock energy [T2l l3T] . 

Figure [8a| demonstrates on the example of the SbH molecule the improvement of the IPEA 
(version A) ground state success probability during the ASP procedure. The dependence 
of the energy gap (AE) between the ground and the first excited state on the adiabatic 



transition parameter s is shown in Figure 8b Although AE is getting close to for r = 8.0 gq 
and s going to 1, the ground state is becoming degenerate at this internuclear distance and 
this fact thus does not infiuence the IPEA success probabihty. We will discuss the relativistic 
qFCI method and its application to the SbH spin-orbit sphtting in more detail in Section [V] 
Recently, the ASP of the hydrogen molecule ground state has been realized experimentally 
on a NMR quantum simulator [2T] . 



C. Controlled "time propagation" 

To study the overall scaling of the qFCI algorithm, one must decompose the only multi 
qubit gate from Figure [SJ i.e. controlled powers oi U = e*"^^ to elementary one and two 
qubit gates. 

For this purpose, it is convenient to express the electronic Hamiltonian in the second 
quantized form [35] 



H 



1 ^ 

y^ hpqalttq + ^ ^ gpqrsala\asar = ^ ^x, (17) 



2 

pq pqrs X=l 



where hpq and gpqrs are one- and two-electron integrals and ol and hp are fermionic cre- 
ation and annihilation operators. The whole summation is formally expressed as a sum of 
individual terms hx- The molecular integrals hpq and Qpqrs can be efficiently precalculated 
on a conventional computer [SB] and represent a classical input to the quantum algorithm. 
In non-relativistic case, they are real-valued, while in relativistic case, they are in general 
complex. 

Since the creation and annihilation operators generally do not commute, the exponential 
of a Hamiltonian cannot be written as a product of the exponentials of individual hx-, but a 
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Figure 8: Adiabatic state preparation (ASP) of the SbH ground state for different 
internuclear distances, (a) Dependence of the IPEA success probabihty on time during the 

ASP (1000 hE^^ ^ 10"^'^ s). Sohd hues correspond to the success probabihties, 

KV'ASplV^exact)!^ ' (0.81, 1) interval is coloured. Reprinted from ^T\. (b) Dependence of the 

energy gap between the ground and the first excited state on the adiabatic transition 

parameter s. 



numerical approximation must be used [2]. The first-order Trotter approximation [37] can 
be expressed as 



JtH 



J-^Y.x^i'hx 
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X=l 



/hxr/N 



N 



+ 0{t''/N). 



(18) 



By choosing A^ > {r^/e), we can implement U within an error tolerance of 0{e) using 
C(L(r2/e)) particular terms e^>'xr/N^ 

Before discussing the factorization of these terms to elementary quantum gates in Section 
IIICl, we would like to mention the implementation modification we use [SO]- Two more 
external inputs are necessary in our case. These are maximum (-Emax) and minimum (-Emin) 
energies expected in the studied system and we use U in the form 



U 



JriEn 



-H) 



For r, it holds 



27r 



-En 



E„ 



(19) 



(20) 



and the final energy is obtained according to the formula 

^ ^ 27r0 
R = R - 



(21) 

The modification mentioned above assures (j) to be in the interval (0, 1). 

ii^min and -Emax cau in fact be chosen arbitrarily, but one must be sure that the calculated 
energy is within this interval, otherwise one would end up with a nonphysical energy, due 
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to the periodicity of e^'^^'^. The maximuin energy can be e.g. the upper bound provided by 
any classical variational {polynomially scaling) method, techniques for calculation of lower 
bounds [38n40j can on the other hand give the minimum energy. The smaller the interval 
between them is, the less iterations of IPEA are necessary for desired precision of E. 

Taking U in the form (19) does not pose any difficulties indeed and as the following circuit 



equality shows, just one more one-qubit rotation is needed. 



jT(En 



-H) 



-irH 



1 

e*^^" 



1. Decomposition of unitary propagator to elementary quantum gates 

The decomposition of the unitary propagator e*^^ to elementary quantum gates [H |2U Hi] 
proceeds in the following manner. First, the Jordan- Wigner transform [42j is used to express 
the fermionic second quantized operators in terms of Pauli a matrices. The Jordan-Wigner 
transform has the form 




(22) 



where a± = 1/2 (cr^^ ± icTy) and the superscript denotes the qubit on which the matrix 




operates. The Hamiltonian (17) then can be rewritten using strings of a matrices. Finally, 



the exponentials of these strings are build up from single-qubit gates and controlled NOT 
operations (CNOTs) [25] . 

We will demonstrate this approach on the one-electron part of the Hamiltonian (with 
complex- valued molecular integrals) 



111 / ^ '^pq'^p'^q / ^ '^ppO'pO'p ~r / ^ yl^pqO'pOiq \ llqpOjqO'pj ■ 

pq pp p>q 

Employing the Jordan-Wigner transform, the diagonal terms can be written as 



(23) 



rippCipCLp 



2 ^ 



<T^ 



(24) 



where h^ is the real part of hpp [/i^p (the imaginary part) is equal to zero due to the Hermicity 
of H] . For the exponentials holds 



-,ihxr/N __ ihppalapT/N 



1 

Q gihppT/N 



(P) 



(25) 



The superscript (p) at the matrix denotes the qubit on which the one qubit gate operates. 
Similarly, the off-diagonal terms read 
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+ 



"'pq 



"'pq 



a; 



cr': 



I, i 7 t 

ripqdpClq "T riqp(Xq(Xp 






(j:. 



+ 



(26) 



where (tI'^'^ represents the direct product 



o 



p^q 



a 



P+i 



al+^® 



a 



q-2 



a 



<?-! 



(27) 



Note that (26) contains the four aforementioned strings of a matrices. 

The exponential of the string of az matrices exp[ir(cr2 ... cr^)] is in fact diagonal in 
the computational basis with the phase shift e^*"^ on the diagonal. The sign of this phase 
shift depends on the parity of the corresponding basis state ( "+" if the number of ones in 
the binary representation is even, "-" otherwise). The exponential can be implemented with 
the following circuit [25j 



-^ 



-^ 



-^ 



(28) 



-^ 



where for the z-rotations holds 



^ Rz{-2t) -^ 



RM 



^-ie/2 
e^^/2 



(29) 



and CNOTs assure the correct sign of the phase shift according to the parity of the state. 
Due to the following change-of-basis identities [25j 



ay = YazY\ 



(30) 

(31) 



where 



the exponentials 



Y = R^{-it/2) 



1 i 
V2\^ 1 



(32) 



exp 



exp 






pq 

2N 



< ® {ar') ® <yl 
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exp 



exp 



ihLr 



PQ 



(J: 



L 2N 

2N 



I ® [aV) ® al 



■o:. 



{oT") ® ^: 



can be implemented with the following circuit pattern, 



V 

p — 1 



^t 



-^ 



-^ 



A 



(33) 



(34) 



g + 1 



? -5t 



^ Rz{Q) -^ 



B 



where A and B are for the individual exponentials (33) equal to {H,H}, {Y,Y}, {Y,H}, 
and {H,Y}, respectively, and 9 to —h^gT/N, —hp^r/N, —hp^r/N, and hp^r/N, respectively. 
Note that although the two strings of a matrices in the first parenthesis in (26) commute as 



do the two strings in the second parenthesis, they do not commute mutually. This, however. 



is not a complication since the Trotter approximation ( 18 ) must be employed anyway. 



We have demonstrated the decomposition technique for the direct mapping approach on 
the one-electron part of the Hamiltonian. The procedure for the two-electron part is more 
elaborate, but completely analogous and we refer an interested reader to [21], where all the 
cases are presented in a systematic way. 

The overall scaling of the algorithm is given by the scaling of a single controlled action 



of the unitary propagator without repetitions enforced by the Trotter approximation ( 18 ) 



These repetitions increase only the prefactor to the polynomial scaling, not the scaling itself. 
Also the required precision is limited, about 20 binary digits of (p are sufficient to achieve 
the chemical accuracy [30] . 

The single controlled action of the exponential of a one-body Hamiltonian (|23| results in 



0{n ) scaling: there are 0{n ) different hpq terms and each of them requires 0{n) elementary 



quantum gates [see the circuit (34)]. Since the same decomposition technique apphed to the 
two-electron part of the Hamiltonian gives rise to similar circuit patterns p4], each term 
gpqrs requires 0{n) elementary quantum gates as well (this in fact holds for general m-body 
Hamiltonians [S]). The total scaling is thus 0{n^) [20l |2^ . where n is the number of 
molecular spin orbitals (or bispinors in the relativistic case [31]) and the qFCI achieves an 
exponential speedup over the conventional FCI. This speedup is demonstrated in Figure [9j 
At this point, we would like to make few remarks. Firstly, we assumed that the initial 
state preparation is an efficient step, as was already mentioned. Secondly, when a quantum 
chemical method with a scaling worse than 0{rv') is used for calculation of the initial guess 
state on a conventional computer, then this classical step becomes a rate determining one. 
Besides this, the classical computation of the integrals in the molecular orbital basis scales 
as 0{n^) (due to the integral transformation) as well. We also assumed noise-free qubits 
and thus did not take into account any quantum error correction [13]. Clark et al. stud- 
ied the resource requirements for a similar, but fault-tolerant computation of the ground 
state of a one dimensional transverse Ising model [H] on a proposed scalable quantum com- 
puting architecture [15]. They showed that due to the exponential scaling of the resource 
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Figure 9: The exponential speedup of the qFCI over the FCI. In case of the FCI (blue), 

dependence of the number of Slater determinants in the FCI expansion on the number of 

basis functions is shown. In case of the qFCI (red), dependence of the number of one and 

two qubit gates needed for a single controlled action of the unitary operator on the number 

of basis functions is presented. Points in the graph correspond to the depicted molecules 

(hydrogen, methylene, methane, ethane, and benzene) in the cc-pVDZ basis set. Reprinted 

from 1301. 



requirements with the desired energy precision as well as due to the Trotter approximation 
employed, an elaborate error correction is required, which leads to a huge increase of a 
computational time. They also gave the values of the experimental parameters (e.g. the 
physical gate time) needed for acceptable computational times. However, the question of 
reducing the resource requirements needed for fault tolerant qFCI computations is still open 
and undergoes an active research. 

IV. APPLICATION TO NON-RELATIVISTIC MOLECULAR HAMILTONIANS 



We will demonstrate an application of the qFCI method to non-relativistic ground and 
excited state energy calculations on the example of the methylene molecule j30j . 



A. Example of CH2 molecule 



Methylene molecule (CII2) in a minimal basis set (ST0-3G) is a simple, yet computa- 
tionally interesting system suitable for simulations and performance testing of the qFCI 
method. CH2 is well known for the multireference character of its lowest-lying singlet elec- 
tronic state (a ^Ai) and is often used as a benchmark system for testing of newly developed 
computational methods (see e.g. [iGl - liQ] ). In the ST0-3G basis set, the total number of 
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molecular (spin) orbit als is 7(14) which means 15 qubits in the direct mapping approach (one 
qubit is needed in the read-out part of the register). Since classical simulations of 15-qubit 
computations are feasible, the CH2 molecule serves as an excellent candidate for the first 
benchmark simulations. 

We simulated the qFCI energy calculations of the four lowest-lying electronic states of 

were studied: C-H 
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CH2: X ^Bi, a ^Ai, b ^Bi, and c ^Ai. Two processes shown in Figure 

bond stretching (both C-H bonds were stretched, Figure 10a), and H-C-H angle bending 



for a Ai state (Figure 10b). These processes were chosen designedly because description of 



bond breaking is a hard task for many of computational methods and H-C-H angle bending 
since the a ^Ai state exhibits a very strong multireference character at linear geometries. 

Our work followed up the work by Wang et al. [T3] who studied the influence of initial 
guesses on the performance of the quantum FCI method on two singlet states of water 
molecule across the bond-dissociation regime. They found out that the Hartree-Fock (HF) 
initial guess is not sufficient for bond dissociation and suggested the use of CASSCF method. 
Few configuration state functions added to the initial guess in fact improved the success 
probability dramatically. 
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Figure 10: (a) Energies of the four simulated states of CH2 for the C-H bond stretching, ro 
denotes the equilibrium bond distance, (b) Energy of a ^Ai state of CH2 for the H-C-H 
angle bending, a denotes the H-C-H angle. Figures reprinted from 



We therefore also tested different initial guesses for the qFCI calculations. Those denoted 
as HF guess were composed only from spin-adapted configurations which quajitatively de- 
scribe certain state. Here, for simplicity, we will present the results just for the X ^Bi ground 
state described by the configuration (lai)^(2ai)^(162)^(3ai)(16i) and a ^Ai state character- 
ized by (lai)^(2ai)^(162)^(3ai)^ configuration. Initial guesses denoted as CAS{x,y) were 
based on complete active space configuration interaction (CASCI) calculations with small 
active spaces, which contained x electrons in y orbitals. Initial guesses were constructed only 
from the configurations whose absolute values of amplitudes were higher than 0.1. Those 
constructed from the configurations whose absolute values of amplitudes were higher than 
0.2 are denoted as CAS{x,y), tresh. 0.2 guess. All the initial guesses were normalized before 
the simulations. 
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In our simulations, similarly as in [12], the exponential of a Hamiltonian operator was 
implemented as an n-qubit gate. This is motivated by the fact that once the decoherence 
is not involved in the model, the exponential of a Hamiltonian can be obtained with an 



arbitrary precision simply by a proper number of repetitions in (18). But we would like to 
note that Whitfield et al. [21] among others also numerically studied the length of a Trotter 
time step needed for a required energy precision on the example of the Helium atom. We 
ran 20 iterations of the IPEA which correspond to the final energy precision ^ 10~^ Eh- 
All the other computational details including the exact definition of the CA spaces can be 
found in [301. 
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Figure 11: Success probabilities of the A version of IPEA for the a ^Ai state with HF and 
CAS (2,2) initial guesses, tresh 0.2 means that only configurations with absolute values of 
amplitudes higher than 0.2 were involved in the initial guess, a denotes the H-C-H angle. 

Reprinted from [30] . 



Figure [TT] presents the results of the A version of IPEA for the angle bending process. 
The overlap and scaled overlap of the initial HF guess wave function and the exact FCI wave 
function is shown as well as the success probabilities for the HF and CAS (2,2) initial guesses 
and dotted line bounding the safe region. The results numerically confirm that the success 
probabilities always lie in the interval | (V'imt|V'exact)| ■ (0.81, l), depending on the value of the 



remainder S (12). This algorithm can be safely used when the resulting success probability 



is higher than 0.5 (as it can then be amplified by repeating the whole process). When going 
to the linear geometry, where the a ^Ai state exhibits a very strong multireference character 
and the restricted HF (RHF) description is no more adequate, the CAS (2,2) initial guess 
improves the success probabilities dramatically. Moreover, these initial states correspond to 
only two configurations and are thus very easy to prepare [3]. 

Performance of the B version of IPEA with the HF and CAS (2,2) initial guesses is 
illustrated in Figure 12 As can be seen, in the situations where the particular initial guess 
state can be used, few repetitions are enough to amplify the success probability to unity. 
The HF initial guess is again not sufficient for the linear and close to linear geometries. 

Results corresponding to the a ^Ai and X ^Bi states for the C-H bond stretching are 
summarized in Figures [13] and [14] Figure [13] presents the performance of the A version 
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Figure 12: Success probabilities of the B version of IPEA with (a) HF and (b) CAS (2,2), 

tresh. 0.2 initial guesses and different number of repetitions of individual bit measurements 

for the a ^Ai state, a denotes the H-C-H angle. Figures reprinted from [30] . 
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Figure 13: Success probabilities of the A version of IPEA for (a) a ^Ai state, (b) X ^Bi 
state of CH2. Different initial guesses were used, tresh 0.2 means that only configurations 
with absolute values of amplitudes higher than 0.2 were involved in the initial guess, tq 
denotes the equilibrium bond distance. Figures reprinted from 



of IPEA. When going to more stretched C-H bonds, the RHF initial guess again fails. 
The CAS (2,2) initial guess improves the success probabilities in case of a ^Ai state near 
the equilibrium geometry but in the region of more stretched C-H bonds it also fails. In 
this region, CAS(4,4) initial guess must be_used [CAS(4,4), tresh. 0.2 guess is sufficient]. 
The situation is even more difficult for the X ^Bi state when the C-H bonds are stretched. 
Here, even the CAS (4,4) initial guess fails and bigger active space - CAS (4,5) - must be used 
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Figure 14: Success probabilities of the B version of IPEA with "best" initial guesses and 
different number of repetitions of individual bit measurements, tq denotes the equilibrium 
bond distance, (a) a Mi state, CAS(4,4), tresh. 0.2 guess, (b) X ^Bi state, CAS(4,5), 
tresh. 0.2 guess. Figures reprinted from 



for initial guess state calculations. However, apart from the CAS size, initial guess states 
always contained at most 12 configurations, but usually 8 or even less for nearly dissociated 
molecule. This observation is in agreement with the results of [13], where few configuration 
state functions added to the initial guess improved the success probability dramatically. 

Figure 14 presents the performance of the B version of IPEA for the "best" initial guesses 
in terms of price/performance ratio. The evident result is again that relatively small number 
of repetitions (~ 31) amplifies the success probability nearly to unity. This is very important, 
because the B version of IPEA seems to be a better candidate for the first real larger-scale 
qFCI calculations, since it does not require a long coherence time. 



V. EXTENSION TO RELATIVISTIC MOLECULAR HAMILTONIANS 



So far, we concerned non-relativistic computations only. But it is well known that rela- 
tivistic effects can be very important in chemistry. In fact, accurate description of systems 
with heavy elements requires adequate treatment of relativistic effects [50] . The most rigor- 
ous approach [besides the quantum electrodynamics (QED) which is not feasible for quantum 
chemical purposes] is the four component (4c) formalism [JT] . Very recently, we therefore de- 
veloped the qFCI method for 4c relativistic computations [HT] and the details of this method 
are the subject of this section. 

The 4c formalism brings three major computational difficulties: (1) working with 4c 
orbitals (bispinors), (2) complex algebra when molecular symmetry is low, and (3) rather 
large Hamiltonian matrix eigenvalue problems [due to larger mixing of states than in the 
non-relativistic case] . All of them can nevertheless be solved on a quantum computer. Before 
discussing how this is done, we would like to note that in our work, we restricted ourselves 
to the 4c Dirac-Coulomb Hamiltonian (with the rest mass term mc^ subtracted): 
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(3' = l3-h, (3 
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where ak [k = x, y, z) are Pauli matrices and I2 the 2x2 unit matrix. This type of 
Hamiltonian covers the major part of the spin-orbit interaction (including two-electron spin- 
own orbit) and also scalar relativistic effects. It is in fact without loss of generality sufficient 
for our purposes since going to Dirac-Coulomb-Breit Hamiltonian |5l] correct to 0{c~'^) 
is conceptually straightforward as the inclusion of the corresponding integrals requires a 
classically polynomial effort. 

We will start the description of the algorithm with a mapping of a relativistic quantum 
chemical wave function onto a quantum register. We have already briefly mentioned this 



topic in Section III We in fact do not need any qubits to represent positronic bispinors in 



the direct mapping due to the no-pair approximation [51] , which is widely used in relativistic 
quantum chemistry. In this approach, one actually builds up an A^-electron wave function 
only from Slater determinants containing positive-energy bispinors. This procedure neglects 
all QED effects, but it is justifiable at the energy scale relevant to chemistry. Moreover, 
because of the time-reversal symmetry of the Dirac equation, bispinors occur in degenerate 
Kramers pairs ^T\ denoted A and B (analogy to a and (3 spin in non-relativistic treatment). 
The direct mapping thus employs one qubit for bispinor A and one for B. The 4c character 
of molecular bispinors does not complicate the approach substantially, since the Hartree- 
Fock (HF) part of a calculation is done on a classical computer and only the exponentially 
scaling FCI part on a quantum one. 

Assuming the no-pair approximation and the empty Dirac picture, the relativistic Hamil- 
tonian has the same second quantized structure (17) as the non-relativistic one. The only 
difference is that the molecular integrals hpq and gpqrs are now in general complex (and 
have thus lower index permutation symmetry). This is actually no difficulty for a quantum 
computer, since our working environment is a complex vector space of qubits anyway and 
we do the exponential of a complex matrix even if the Hamiltonian is real. Moreover, on 
the example of the one-electron part of the Hamiltonian [see (26) and the circuit (34)], one 
can see that the decomposition of the unitary propagator e*"^^ with complex molecular inte- 
grals requires twice as many gates compared to real- valued Hamiltonians pi], while complex 



arithmetic on a classical computer requires four times more operations. 

The last of the aforementioned difficulties of the 4c formalism is the size of a Hamiltonian 
matrix eigenvalue problem. When we put the double group symmetry aside and employ 
Kramers restricted (KR) approach, the relativistic Hamiltonian, unlike the non-relativistic 
one, mixes determinants with different values of the pseudo-quantum number Mx 



Mk = 1/2{Na - Nb) 



(36) 



where Na and Nb denote the number of electrons in A and B Kramers pair bispinors (in 
the non-relativistic case, Mk corresponds to Ms and Hamiltonian is block diagonal in Ms). 
Therefore, in the case of a closed shell system with n electrons in m molecular orbitals 
(bispinors), the size of the relativistic Hamiltonian matrix (number of Slater determinants 
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in the FCI expansion) reads 



x=0 



'v.=..-E(: „'" -rr). (37) 



The last equahty is due to the Vandernionde's convolution [52j. When compared to the 
non-relativistic case ([I]) and using the Stirling's approximation in the form 

1 
In m! ~ -In {27im) + mln m — m for m — t- oo, (38) 

the ratio between the relativistic and non-relativistic Hamiltonian matrix sizes is given by 
the expression 



where we set m = fc ■ n. 

On a quantum computer on the other hand, this problem does not occur. When employing 
the direct mapping which maps the whole Fock space of the system on the Hilbert space of 



qubits, the Hamiltonian (17) in fact implicitly works with all possible values of M^- The 



scaling of the relativistic qFCI method is therefore the same as the non-relativistic one, 



namely 0{m^), where m is the number of molecular orbitals/bispinors (see Section IIIC 1). 
It is quite remarkable and deserves an emphasis that the relativistic qFCI method not 
only achieves an exponential speedup over its classical counterpart, but, as we have just 
discussed, also has the same cost (in terms of scaling) as its non-relativistic analogue. 



A. Example of SbH molecule 

We will demonstrate the performance of the relativistic qFCI method on the example 
of the SbH molecule. This molecule is of particular interest to us, since its non-relativistic 
ground state ^S~ splits due to spin-orbit effects into X O"*" and A 1 states. In the approximate 
Aw-projection, these states are dominated by o"i/2'''"i/2^3/2 ^^^ '^\i2^\i2^\i2 configurations. 
The splitting is truly of "molecular nature" as it disappears for dissociated atoms and its 
experimental value is Ai^so = 654.97 cm~^ [SH] . 

For this reason, we simulated the SbH bond dissociation process. Simulated potential 



energy curves of both states are shown in Figure 15, Since we employed rather large basis 
sets of triple-C quality, we could not manage to simulate the FCI calculations with all 
electrons. We instead simulated general active space (GAS) KRCI computations [51] with 
the occupation constraints that give rise to CI spaces of approximately 30000 determinants. 
Definition of the GA spaces together with all the computational details can be found in [31] . 

We worked solely with a compact mapping employing the double-group symmetry {C\^ 
and exponential of a Hamiltonian was again simulated as an n-qubit gate. We used the 
DIRAC program [SS] for calculations of Hamiltonian matrices and ran 17 iterations of the 
IPEA which correspond to the final energy precision ?^3.81 x 10"^ E^. 

Based on our KRCI setup we obtain a vertical Aii^so of 617 cm~^. Success probabilities 
of the algorithm with the HF initial guesses (represented by the dominant configurations of 



both states) are presented in Figure 16 They correspond to the A version of IPEA. 
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Figure 15: Simulated potential energy curves of ground (0"*") and excited (1) states of SbH, 
and spin-orbit energy splitting. Reprinted from [3T] . 
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Figure 16: SbH (a) ground (0+) and (b) excited (1) state qFCI (IPEA version A) success 
probabilities corresponding to the HF initial guesses. 



Ground state success probabilities confirm that relativistic states have, due to near de- 
generacies caused by the spin-orbit coupling, often a stronger multireference character than 
non-relativistic ones. The upper bound of the success probability is less than 0.7 even for 
the equilibrium geometry and the lower bound is higher than 0.5 only up to 4.8 Oq. The 
success probabilities of the A 1 state are higher and the HF initial guesses can be in a 
noise- free environment used up to 6 ao- For longer distances, initial guesses from more 
accurate polynomially scahng relativistic methods must be used (analogously to the non- 
relativistic example, e.g. the relativstic 4c CASSCF method with a small orbital CAS). 
The low ground state success probabilities can be also improved by the ASP procedure as 
is shown in Figure |8l 
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VI. CONCLUSIONS 

In this article, we have reviewed quantum algorithms for the FCI energy calculations. 
Starting with the seminal paper by Aspuru-Guzik et al. [12j and ending with our most recent 
work on relativistic generalization of the algorithm [31], we gave a detailed description of the 
qFCI method which is applicable to the ground as well as excited state energy calculations 
in the non-relativistic and also relativistic regimes. Both variants share the scaling in the 
number of molecular orbitals/bispinors m, namely 0{m^) and exhibit thus an exponential 
speedup over their classical counterparts. We have demonstrated their performance by 
numerical simulations of the CH2 and SbH energy calculations. The results indicate that 
the ground and excited state energies at equilibrium geometries are accessible with HF initial 
guesses, which are easy to prepare. On the example of the CH2 molecule, we have shown 
that CASCI initial guess states with small complete active spaces composed of relatively 
few configurations (~ 10) are sufficient even for a nearly dissociated molecule to achieve 
the probability amplification regime of the IPEA algorithm. The first proof-of-principle 
experimental realizations recently achieved by several groups [2nH22] are very promising and 
confirm the real applicability of the method. 
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